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TECHNICAL MEMORANDUM X-53951 


MIE SCATTERING: A COMPUTER PROGRAM AND AN ATLAS 

SUMMARY 


A computer program for calculating the Mie Scattering parameters is 
presented with applications in distributions of particle size and effects of 
different lens apertures of a detector. The computer program, * MIESCA, 
calculates the amplitude functions, efficiency factors, intensities, and polari- 
zations. The program was adapted for magnetic tape filing in the Univac 1108 
and for plotting intensities and polarization versus scattering angle on the 
SC 4020 plotter. An Appendix and an Atlas of scattering intensity graphs are 
included. This report provides a basis for a program to evaluate character- 
istics of the radiation scattered by a dilute cloud of particles. 


INTRODUCTION 


Electromagnetic radiation passing through a dilute cloud of particles is 
scattered and absorbed by the particles in the cloud. The characteristics of 
the scattered radiation are determined by the wavelength, X, of the incident 
radiation, the complex refractive index, m, of the medium, and the size and 
shape of the discrete particles in the medium. Therefore, information about 
the state of the medium can be provided by the measurement and proper 
interpretation of the characteristics of the scattered radiation. The numerical 
determination of the characteristics of the scattered radiation for given particle 
models is of prime importance. This report provides the first step in a 
program for the evaluation of characteristics of the radiation scattered by a 
dilute cloud of particles. Programs for scattering by a single spherical 
particle (Mie Scattering) are presented. 

T-027 (ATM Contamination Measurement) and S-073 (Gegenshein/ 
Zodiacal Light) are two photometric experiments assigned to the first AAP 
flights. The T-027 objective is to study the spacecraft contaminant cloud by 
observing the scattered light of the particles. The S-073 objective is to study 
the interplanetary dust by scattered light. The scattering programs presented 
and to be developed will be available for use in data analysis for the T-027 and 
S-073 flight experiments. 


* The basic program was developed by Brown Engineering Company under 
NASA Contract NAS8-20166. 



MIE SCATTERING 


Theory 


In the T-027 and S-073 AAP experiments, light scattering by particles 
will prove to be important in data analysis. Therefore, this initial report 
provides an Atlas of Mie Scattering diagrams for different particle sizes to 
aid in the identification of possible contamination for T-027 and in synthesizing 
a zodiacal light model for S-073. Light scattering by spheres of arbitrary size 
and refractive index are described by Mie Scattering. The angular distribution 
of intensity of scattered radiation was obtained by Mie, and independently by 
others, who successfully applied Maxwell's equations with the appropriate 
boundary conditions. The basic expressions for the radiation scattered by a 
sphere of radius a, and of a material with a complex refractive index m, 
are [ij: 
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where P n (x) is a Legendre polynomial of order n, and x = -cos (0) . A^ 

and B^ are amplitude functions (linear combinations of Ricatti-Bessel, 

Ricatti-Hankel functions and their derivatives, with arguments a = (27ra)/A 
and ma). The physical significance of ij and i 2 , the intensity functions, is 
that they are the intensities of the two incoherent plane-polarized components 
of light scattered by a single particle (Fig. i) . When this particle is illuminated 
by randomly polarized light of unit intensity, then (A 2 ^)/ (87T 2 ) and (A 2 i 2 ) / (87T 2 ) 
are the radiant intensities of the scattered components with electric vectors 
perpendicular and parallel, respectively, to the plane of observation. The 
differential scattering cross section per particle is 
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Figure 1. Physical representation of the intensities. 


di + y 

Q is the scattering coefficient or efficiency, defined by the ratio of the 

SCcl 

scattering cross section to the geometric cross section of the particles. 

(The scattered wave intensity is obtained by integration of the radial component 
of the Poynting vector of the scattered wave over the surface of a sphere 
surrounding the particle. ) The angular distribution function, f (6 ) , the 
fraction of scattered radiation directed into a solid angle in the direction of 6 , 
is defined as 



Substituting the expression for f(0) into the expression for the differential 
scattering cross section yields: 
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The Computer Program 

The form of the general formulas is not amenable to computer computa- 
tion, so numerical methods must be employed [i]. The most fundamental 
factors in the Mie solution are the amplitude functions, which can be calculated 
by the following formulas: 
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where 

ib is a Ricatti-Bessel function, of order n. and 
n 

£ is a Ricatti-Hankel function, of order n. 
n 

Also, 

a = — — and /3 = ma = a(m*-ik) , 

A 

where the basic parameters are defined: 

a = radius of a spherical particle 
A = wavelength of incident radiation 
m = complex refractive index 
m* = real part of m 
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k = imaginary part of m (extinction coefficient of the particle 
material) 

i = (-0* 

From the amplitude functions, the efficiency factors for total extinction, 
scattering and absorption are computed: 
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The computer input is composed of the scattering angle, the refractive 
index of the medium, and a . The program computes and prints out the 
efficiency factors, both intensities ij and i 2 , their sum, and the polarization. 
(The polarization is defined as: i* - i 2 /ii + ^ ) An example of the computer 
output is given in Table 1, for which a = ALT = 5. 0, 0 = XTH1 = 0 degree 
(5 degrees) 180 degrees, and m* = EM = 1.33. An SC 4020 plot for these 
values is also included in the Atlas. Appendix A gives the details for the 
graphical display using the Univac 1108 — SC 4020 plotter which was used to 
produce the Atlas (Appendix C). Appendix A gives the 1108 user information 
concerning the Mie computer programs. 


APPLICATIONS 


Dilute clouds of particles are normally polydispersed and are observed 
through finite apertures; therefore, integrations over 0 and a are necessary 
for data analysis. This section considers integrations over particle size 
distribution and finite angular apertures. 
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TABLE i. SAMPLE OF MIESCA PRINTOUT FOR a = 5. 0, m = 1.33 
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Distributions 


The zodiacal cloud is composed of many different sizes of particles, 
whose distribution is generally assumed to be of an exponential form: 



where a and a are as previously defined, and p is dependent on a^, 

a , and some "breakpoint" a_ between them [2] . Other types of distri- 
max B 

buttons, such as aerosol distributions, will be considered for T-027. 


The first distribution examined was the exponentially decreasing 

-O' -oi 

function e . The choice of e is convenient because, 


/ e~ a da = 1 
0 

So, the percentage of each size particle in the distribution can be found easily 
on the basis of its fractional part of the whole. 

When calculating an intensity for particles with a size distribution, the 
different intensities for each a must be weighted according to the distribution. 
For a particular 0 , 
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Then, this must be normalized to give 
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By Gauss-Legendre Quadrature, this can be numerically integrated in a 

Fortran program. (See Appendix B. ) An example for m = i. 3 - 0. Oi, 

a . = 0. 1 and a = 10. 0 is in Table 2. 
mm max 

Apertures 

The angular aperture of the objective lens of a photometric system 
distorts the observed scattering distribution. This effect is being studied in 
connection with the T-027 data analysis. It is impossible to record the intensity 
at a single 0, just as it is impossible to measure the integral or "area" under 
a point. An integral must be taken over some interval. The interval in this 
study is the lens aperture. The intensity is calculated over the whole aperture 
and the "integrated" intensity distribution is determined. Effects of changing 
the angular aperture of the flux collecting device must be studied to insure 
correct interpretation of measured scattering distributions. Widening the 
lens aperture smoothes the intensity distributions and reduces the number of 
maxima. 

The flux, $, is the integral of the intensity over the solid angle. 

Gucker and Tuma [3] have developed a method of modifying the flux equation, 
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for computation, where 17 = angle of observation, co is the angle of aperture, 
and 0 is the scattering angle (Fig. 2). After normalizing (by dividing 

/ A \ 2 

through by ) > the flux equation becomes: 
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TABLE 2. NORMALIZED DISTRIBUTION FOR m = 1. 3 - 0. 0i 
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where 
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The angular apertures are divided into n strips, which are, in turn, 
subdivided into b substrips. Table 3 was produced from the computer 
program using an aperture of 2co = 6 degrees, an observation angle rj, 
n = 5, b = 4, and with the Mie program (a = 10. 0, m = 1. 33) . Figure 3 
shows implicitly the aperture smoothing effect for these values (cf. Atlas 
for a = 10) . 

More applications for the Mie program will be developed in the future. 
These will include a further analysis of the smoothing effects of both the 
polydispersions and the aperture. Other computer programs to be developed 
will be concerned with 

1 . Zodiacal models developed in full 

2. Contaminant cloud models developed in detail to predict distributions 
and intensities 

3. Problems associated with nonspherical particles. 
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TOTAL INTENSITY 



150 ° 155 ° 

SCATTERING ANGLE (DEGREES) 


Figure 3. Intensity comparison of a 6-degree aperture with 
that of a 0-degree aperture. 
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APPENDIX A 


M IE COMPUTER PROGRAMS 
File Tape on Univac 1108 

A more efficient use of computer time has been realized by filing large 
subroutines on tape for use in the Univac 1108 computer. Then, the calling 
program (e. g. , an integration) , can have a minimum execution time with no 
computer printout. The resultant control cards which are used for filing the 
Mie program on magnetic tape are given in Program List 1 where MIESCA 
is the Run Identification, MIE is the filed subroutine, MIETAP is the name of 
the newly formed tape, 000000 is the identification (ID) number, and 5 is the 
maximum number of minutes of computer time. 
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PROGRAM LIST 1. FILE SUBROUTINE MIE ON MIETAP 


3 RUN » A MI£SCA.0ODeeO*NADrNtBIN219*5 
SASG,T MIETAP,T, 32074 

9MSG,W PLEASE PUT WRITE RING IN TAPE 17074 
3F0R,IS HIE 

SUBROUTINE MIE t XTH 1 , EM , AL T , X li *X 12 » X12 tPOL AR ) 

DIMENSION EM 12), BETC2) ,XJNB<7) , X JN8P 1 2 ) * AUS3 (2 I , RB SB C21.RBSBPC2J 
1 ,RHNK< 2) ,RHNKPf2),Pl C 2 ) , P2 C 7 I ♦ P3 ( 2 ) t XNUMA i 2) tXNUMB (2 ) ,XDENA(2> , 

1 XDENR < 2) » AN ( 2) ,BN C 7 ) ,XTH<37) *PN < 37 > , PNP < 37 ) , SI NR C 3 7) , SINK 37), 
1S2NRI37) ,S2NI(37),S1RC 37>,S2RC 37 > , SI I ( 37 > » S2I( 37), XX <3 7) , 

IAAUKSC 37),AAUX (37) 

COMMON N , XN 9 ALP, BET , AUX1 , AUX 2, AUX3 1 AUX4, AUX5 (2 ) , AUX6 C 2 ) , AUX, 

I AUKS , VX 
1=1 

ALP=ALT 
XTH II) rXTHl 
PI=3 * 1 41 5977 
DEGzPI/l 80 ® 

52 CALL 7XRCFM,ALP,BET»1> 

13 M=1 

XTHCI)=XTHII)*DEG 
56 SI R f I) =Q« 

SI I ( T) =0 » 

S2R 1 1 ) =0* 

11 S?I<I>=0* 

N=1 

AEXTrO* 

ASCATrO* 
l? XNrFLOATIN) 

CALL RECJf X JN, XJNP»X JNI, XJN1P, XJNB ,X JNBP*AUS1» AUS2 ,AUS3) 

AUX1=XJN 
AUX2=XJNI 
A UX3=AUS 1 
AUX4=AUS7 

CALL 7EQUAfAUX5,XJNB? 

CALL 7EQUAC AUX6*AUS3> 

RBS= ALP ♦X JN 

RBSP=XJN+ALP*XJNP 

CALL 7MULT ( BET ,XJNB« RBSB ) 

CALL 7MULTCBET *XJNBP*P1) 

CALL 2ADDIXJNB*P1*RBSBP) 

RHNK ( 1 ) =RBS 

RHNK<2)=f-l. )**N*XJN1*ALP 
RHNKP< 1 ) rRBSP 

RHNKP ( ?)=t i -l »)**N*XJN1P*ALP+ «-!l **N*XJN 1 
CA=( CXN+XN+1 •)*!-! ) **N 1/ I XN* XN*XN ) 

CBr-CA 

CALL ZXR«RBSBP,RBS,PI*1> 

CALL 7XR(RBSB,RBSP,P2* 1) 

CALL 7MULTCEM,P2.P3> 

CALL 7SUB(P1 *P3» XNUMA ) 

CALL 7MULTf£W,Pl»P3) 

CALL 7 SUB I P3 »P2, XNUMB) 

CALL 7MULT(RHNK,RBSBPpPU 
CALL 7MULT(RBSB, RHNKP, P2J 
CALL 7MULT f EM»P2 ,P3) 

CALL 7SUBI PI ,P3,XDENA) 

CALL 7MULTfEM,Pl ,P3) 

CALL 7SUB«P3,P2,XDENB) 

CALL 7XICXNUMA,CA,PI»1) 

CALL 7DIV( Pi »XDENA » AN) 
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PROGRAM LIST i. (Continued) 


CALL 7X1 CXNUMB.CB.P1 »1 I 
CALL 7DIVCP1 .XDENB.BN) 

17 CXNrXN* CXN+1.)*(— 1.1**N 
AEXN-CXN*CANC21— BNC21) 

AEXT-AEXT.-AEXN 

ABSANr ZMOO t AN ) 

ABSBNrZMOD C BN) 

CKN~CXN*XN*CXN+1.)*CXN*1.) 1/ CXN.-XN-.1 .) 

ASCNrCKN*C ABSAN*ABSAN-,ABSBN*ABSBN) 

ASCAT=ASCAT*ASCN 

IF ( ABSAN-. OOOOOO I) 21 .21.23 

21 IFCAB55BN-.OOOOOOU22.22.23 

22 IFCXN-ALP123.23.25 

23 CONTINUE 

18 CONTINUE 

VX=— COSCXTHCI) 1 
IFCN-1 141.41 .42 
42 AUKSSAAUKSCI1 
AUXrAAUXCI) 

41 CALL RECP C XL NP .XLNPP , AUS 1 
PNCIlrXLNP 
PNP C II r XLNPP 
AAUKSCI)rPN(I) 

AAUXCI)=AUS 

XXCI)=XN*CXN+1.)*PNCI1-VX*PNPCI1 
SlNRCIlrANCl )*PNPCI)+BNC 1)*XX< I) 
S1NICI)=ANC2)*PNPCI)+BNC?)*XXCI1 
S2NR CT irBN C 1 >*PNP C II ««N(1) *XXC II 
S2NICI)=BNC2)*PNPCI)+ANC2»*XXm 
SlRCIlrSIRCIH-SlNRCI) 

SlICIlzSlICIH-SlNICI) 

S2R C 1 1 =S?R < 1 1 *S2NR C 1 1 

100 S2ICI)=S2ICI)«-S2NI<I) 

101 N=N*1 

GO TO 12 

25 ASCAT=ASCAT*2./C ALP*ALP) 

AEXT=AEXT*2./C ALP* ALP) 

AABS-AEXT— ASCAT 

Xll=SlRCI)4SlRtIH-SlICI)*SlICI 1 
XI2=S7RCI)*S2RCI)*S2ICI)*S2ICI 1 
116 POLAR-CXIl— XI21/CXI1+XI2) 

X12=XI1*XI2 

RETURN 

END 

3F0R.IS RECP 

SUBROUTINE RECP(PN.PNP.AUS) 

DIMENSION BET C 2) 

COMMON N.XN. ALP. BET, AUX1.AUX2. AUX3 . AUX4. AUX5 C2 )♦ AUX6 C2 1 » AUX, 
1AUKS.VX 
IFCN-1 >1 ,1 ,2 

1 PN=VX 
PNPr). 

AUS- CC2.*XN+1.)*PN*VX— XNl/CXN*!. I 
GO TO 3 

2 PNrAlJX 

IF C ABS C VX- 1 . OOOOOO 1—0. 000001 121.21 .22 

22 IF C ABS CVX* 1.000000 1-0. 000001 121,21.23 
71 PNP=CXN*XN,XN)»PN/ CVX+VX) 

GO TO 24 

23 PNPrC— XN*VX*PN+XN*AUKS)/C1 .— VX*VX) 

24 AUS- ( C2.*XN+I.)»VX*PN— XN*AUKS)/CXN>1.) 
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PROGRAM LIST 1. (Continued) 


3 RETURN 
END 

3F0R»IS RECJ 

SUBROUTINE REC J( X JN. X JNP »X JN I, XJN1 P. X JNB »X JNBP » A US It AUS2 * AUS3) 
DIHENSION XJNBC2)* X JNBP C 2) »XJnB(2) ,XJ 1 BC 2 ) t8ET(2). 

1XJ2B(2)»P1 (?)»P2(2)»G1 C2>.02(?>. AUS3 < 2 ) . SIN8 <2 ) . COS8 (2 ) 

COMMON N,XN» ALP.BET.AUX1 tAUX2. AUX3 . AUX4. AUX5 (2 ) » AUX6 C2)t AUX. 
1AUKS.VX 
IFtN-1 >1 »1 7 ? 

1 XJO=SIN< ALP) /ALP 
XJIrXJO/ ALP-COS ( ALP) /ALP 
XJNzXJl 

XJHl:XJO/ALP— XJ1 
XJM2=(-XJMI) /ALP-X JO 
XJNirXJM? 

XJ2=3.»XJ1 /ALP-X JO 
XJNP=(XJO-2.*XJ2J/3. 

XJH3=(-3.)*XJM2/ALP-XJM1 
X JNIPr ( ( -2 • ) *X JH3+X JMl ) / (— 3. ) 

CALL 7SC (BET »SINB* C0S8 » 1 ) 

CALL ZOIV(SINB.BET.XJOB) 

CALL 7DTV(XJ0B»BETrQl) 

CALL ZSC(BET.SINB.C0SB*2> 

CALL 7DIV(C0SBVBET»Q2) 

CALL: ZSUB<01 >Q?<XJ1B) 

CALL ZEOUACXJNBf XJ1B) 

CALL 7XR(XJ1B?3.»P1»I) 

CALL ZOIVIPl »BET »01 ) 

CALL ZSUB ( Q1 »X JOB • X J2B ) 

CALL ZXR(XJ2B.2.tPl»ll 
CALL ZSUB(XJ0B.P1»P2) 

CALL 7XR(P2»3. »XJNBPt2J 

AUS1=XJ2 

AUS2=XJM3 

CALL ZEQUA(AUS3.XJ2B) 

GO TO 3 

2 XJN=AIIX3 
X JN1=AUX4 

AUS1=(2.*XN+1. )*AUX3/ALP-AUX1 

XJNP = (XN*AUX1-(XN+1)*AUS1)/(2.*XN«-1.1 

AUS2=AUX4*(-2.*(XN+1.)+1.)/ALP-AUX2 

XJN1P=(-(XN«-1. >*AUS2 + XN»AUX2>/(-2.*(XN+l .) +1 .) 

CALL ZEQUA ( X JNB « AUX6 ) 

C=2.»XN+1. 

CALL ZXR(AUX6»C»P1*1) 

CALL ZDI V( PI *BET »01 ) 

CALL ZSUB(Q1*AUX5*AUS3> 

CALL ZXR(AUX5tXN.Pl.l) 

CALL ZXR(AUS3.XN+l,*P2.1) 

CALL ZSUB ( Pl't P2* PI ) 

CALL 7XR(P1.C.XJNBP»2) 

3 RETURN 
END 

aFORfIS ZDI V 

SUBROUTINE ZDIV(A.B.C) 

DIMENSION A ( 2) »B(2)»C(2) »BSTAR(2> 

CALL ZCONJ(B.BSTAR) 

CALL ZMULT(A.BSTAR.C) 

C( 1 )=C(l)/(7MOD(B)*ZMOD(B> ) 

C(21rC(2)/(ZM0D(B)*ZM00(B) ) 

RETURN 
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PROGRAM LIST 1. (Continued) 


END 

aroR*is zsur 

SUBROUTINE ZSUBIA.B.CI 
DIMENSION A(2) .0(2>.Cf2> 

Cf l>=ACl hBU) 

C(21rA<2)-B«2) 

RETURN 

END 

aroRsis zaod 

SUBROUTINE ZADO(A.B.C) 
DIMENSION A<2> »B(21*C<2) 
cci)=Acn+BU> 

Cf?5=Al??*B?2) 

RETURN 

END 

aEORtlS ZHUt 

SUBROUTINE ZMULTf A tB *C ) 
DIMENSION A ( 2) »B C2 ) * C ( 2) 

Cf ll=A<l!*Bm~JU2t*BC21 
C < 2) -A < 1 )*B(2>*A(21*B( l) 

RETURN 

END 

3F0R.IS ZCON 

SUBROUTINE ZCON J < A . ASTAR 1 
DIMENSION A ( 2) » AST ARC 2 1 
ASTARmrAei) 

ASTAR«2)r-A«21 

RETURN 

END 

aFOR*is zMon 

FUNCTION ZMOD< A) 

DIMENSION A ( 21 

ZMODrSORTUm*A(l H-A(2)*A<2I) 

RETURN 

END 

9FOR.IS RTHX 

SUBROUTINE RTHXY f RTH t A 1 
DIMENSION RTH (2) .AC21 
AlllrRTHU >*COSCRTHl 21 1 
A C 2 1 -RTH 1 1 1 *SIN I RTH f Z) ) 

RETURN 

END 

3F0R.IS ZXR 

SUBROUTINE 7XR f A *R »C *N ) 
DIMENSION A«21*CI21 
GO TOUt?) »N 

1 CfI)rACll*R 
C ( 2) -A C2 1 *R 
GO TO 3 

2 CClirAUl/R 
C«2irA«2)/R 

3 RETURN 
END 

aFOR»!S ZXI 

SUBROUTINE ZXI C A .XI* Ct N 1 
DIMENSION A < 21 tC C? 1 
GO T0C1.21.N 
2 CU!=-AI?J*XI 
C«?>=AU>*XI 
60 TO 3 

2 Cm=A«21/Xl 
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PROGRAM LIST 1. (Concluded) 


Ct2)=-A(l)/XI 
3 RETURN 
END 

3F0R.IS ZSC 

SUBROUTINE ZSC (A .SIN A, COS A .N J 

DIMENSION A (2) .SINA(7J »COSA<?> ,C1(21 .C2<2> »RTH 1 ( 2) .RTH212) 
RTH1 C! 1=1. 

RTH1 (? ) = A ( 1 1 
RTH211 )=1. 

RTH2 C? )=-A C 1 > 

ER1=EXPC-A(2I» 

ER2=EXP ( A 1 2) ) 

CALL RTHXV(RTHl.Cl) 

CALL RTHXY(RTH2.C2) 

CALL ZXRtClfERl* Cltl 1 
CALL ZXR1C2»ER2»C2»l) 

GO TO ( 1 »2 ) *N 

1 CALL ZSUBJC1 •C2.C1) 

CALL 7XIIC1.2. .SINA.2) 

GO TO 3 

2 CALL 7ADDIC1 »C2.C1> 

CALL ZXR1C1.2. .C0SA.2J 

3 RETURN 
END 

3F0R »IS ZEGU 

SUBROUTINE ZEQUACA.B) 

DIMENSION A ( 2) »B (2 ) 

a 1 1 ) =b ( n 

A ( 2) =B 1 2 J 

RETURN 

END 

aMAP .IS A * R 
IN MIE 

aCOPOUT TPF$. .MIETAP. 

BREMIND MIETAP. 

BERS 

3C0PIN MIETAP. ♦TPFS., 

3REWIND MIETAP. 

3FIN 

3FIN 


For calling MIE out of the tape file [4] , the following sequence is used, 
with tape label 32024, Run ID, and calling program DRYRUN. 

PROGRAM LIST 2. CALL UP FILED MIETAP IN DRYRUN 

@RUN, A DRYRUN, 000000, NAME 
@ASG, T MIETAP, T, 32024 
@COPIN MIETAP. , TPF$. 

@REWIND MIETAP 
@FOR, IS DRYRUN 


Main program 

• * * 
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CALL MIE (XTH1, EM, ALT, XIi,XI2,X12, POLAR) 

e o • 

@MAP, IS A, B 
@X 

IN DRYRUN 
@XQT 
Data 
@FIN 


SC 4020 Plots 

The plots used in the Atlas were done by the SC 4020 plotter called up in 
the Univac 1108 MIESCA program. For the example in the text with a = 5. 0, 
we used Program List 3 for the main program. (Subroutines are the same as 
for Program List 1. ) 

The dimension statements and the data statements supplying the abscissa 
and ordinate labels are labeled A. 


The first card needed for addressing the plotter is the CALL IDENT 
(CAMRAS) , labeled B. Depending on the number in parentheses, the output 
will be on paper or microfilm: 

CAMRAS = 35 film 
CAMRAS = 9 paper 
CAMRAS = 1 film and paper. 

The values for the array to be plotted are generated in the program. The 
CALL SMXYV (MX, MY) is inserted (labeled E) , where MX, MY are fixed- 
point integers which designate whether the linear or nonlinear (logarithmic) 
mode is to be used. 


If [4] MX + 0, MY + 0 
MX + 0, MY = 0 
MX = 0, MY + 0 
MX = 0, MY = 0 


log in X, log in Y 
log in X, linear in Y 
linear in X, log in Y 
linear in X, linear in Y . 


SMXYV must be called each time before there is a change of mode. Now, 
QUKLOG and QUIK3V can be labeled F. The general forms for calling the 
subroutines are: 
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CALL QUIK3V ( L, ISYM, FLDX, FLDY, NP, X, Y) 

CALL QUIK3L (L, XL, XR, YB, YT, ISYM, FLDX, FLDY, NP, X, Y) 
CALL QUKLOG (L, XL, XR, YB, YT, ISYM, FLDX, FLDY, NP, X, Y) 



PROGRAM LIST 3. SC 4020 ADAPTATION OF MIESCA 
(MAIN PROGRAM ONLY) 


SRUN.A MIESCA » 400080 , G ARYHIESC A T 1 » 2 
SFORtlS GARYGA 

C THIS PROGRAM COMPUTES THE MIE SOLUTIONS FOR ARBITRARY PHYSICAL 
C PARAMETERS. 

C 

C DESCRIPTION OF INPUT PARAMETERS - 

C N I TOTAL NUMBFR OF SETS OF INPUT PHYSICAL PARAMETERS 

C EACH SET OF INPUT PHYSICAL PARAMETERS CONTAINS - 

C ALP THE PARAMETER ALPHA 

C EM THE COMPLEX REFRACTIVE INDEX 

C EACH SET OE INPUT PHYSICAL PARAMETERS IS FOLLOWED BY A SET OF 
C OPTION CODFS - 

C N0PT1 ONE DIGIT NUMBER INDICATING HOW THE INTENSITY 

C FUNCTIONS ARE TO BE COMPUTEO 

C 1 - THE INTENSITY FUNCTIONS ARE COMPUTEO FOR 

C SCATTERING ANGLES VARYING FROM 0 TO 180 DEGREES 

C AT 5 DEGREE INTERVALS 

C 2 - THE INTENSITY FUNCTIONS ARE COMPUTEO FOR SPECIFIED 

C SCATTERING ANGLES 

C N0PT2 ONE DIGIT NUMBER INDICATING WHETHER THE AMPLITUDE 

C FUNCTIONS ARE TO BF PRINTEO 

C 1 - YFS 

C 2 - NO 

C N0PT3 ONE DIGIT NUMBER INDICATING WHETHER THE INTENSITY 

C FUNCTIONS ARE TO BE COMPUTEO 

C I - YES 

C 2 - NO 

C NOTE - WHEN N0PT3=2» IT OVERRIDES N0PT1 . N0PT1 MUST BE SET EQUAL 
C TO J IN THIS CASE 

C WHEN N0PT1 =2, ADDITIONAL PARAMETERS MUST FOLLOW THE SAME SET OF 
C INPUT PHYSICAL PARAMETERS - 

C M NUMBER OF SPECIFIED SCATTERING ANGLES 

C NOTE - M MUST BE SMALLER THAN OR EQUAL TO 37 

C XTH A VECTOR OF LENGHTH M WITH XTH(I) EQUAL TO THE I-TH 

C SPECIFIED SCATTERING ANGLE EXPRESSED IN DEGREES 

C 

C DESCRIPTION OF OUTPUT QUANTITIES - 

C AN, BN THE COMPLEX AMPLITUDE FUNCTIONS 

C Xll.Xl? — THE INTENSITY FUNCTIONS 

C AEXT THE EFFICIENCY FACTOR FOR TOTAL EXTINCTION 

C ASCAT THE EFFICIENCY FACTOR FOR SCATTERING 

C AABS THE EFFICIENCY FACTOR FOR ABSORPTION 

C POLAR DEGREE OF POLARIZATION 

C ALL OUTPUT QUANTITIES ARE PRINTFO UNLESS THEY HAVE BEEN RULED OUT 

C BY THE OPTION COOES 

C 

C REMARKS - 

C IN ALL COMPUTATIONS OF INFINITE SERIES, CONVERGENCE IS ASSUMEO 

C TO HAVE BEEN ACHIEVED WHEN BOTH AMPLITUDE FUNCTIONS ARE SMALLER 

C THAN OR EQUAL TO 1.0E-07, AND THE NUMBER OF TERMS COMPUTEO 
C EXCEEDS THE VALUE OF ALP. 

C 

% DIMENSION EM (2 ) , BET ( 2) ,X JNBl 21 »X JNBP <2 I » AUS3 (2 > . RBSB (2 ) » RBSBP< 2) 
1 , RHNK (21 ,RHNKP(2).P1 (2).P2(2).P7(2).XNUMA(2) »XNUMB (2 1 , XDFN A( 2 > , 

1 XDENB 12) »AN{2) ,BN<2> ,XTH(371 ,PN( 77 ) . PNP ( 37 ) . S1NR (37) ,S1NI (37). 
1S2NR (37) . S2NI (37). S1 R( 37). S2R ( 37) , SI I ( 37 ) . S2 1 ( 37 ) , X X(37) , 

1 A AUKS (_37) , A AUX (_37_) 

DIMENSION ARRAY1 (2000) ,BRR A Y 1 ( 2000 ) . CRRA Y 1 (7 DP 0 ) ,DRR A Y 1 ( 70 Of) ) . 

1 ARR A Y7 (7000) , BRR A Y 7 ( 20 00 ) » CRPA Y? ( 7000 ) ,DRR AY?( 7D0D ) 

A DIMENSION AFDX (17) , AFP Y ( 1? ) . BFDY ( 1 ? ) , CFD Y ( 12 ) . OFOY ( 1 2 ) 
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PROGRAM LIST 3. (Continued) 


A 


( OAT A AFDX/72H TMFTA DEGREES 
1 

data afoy/7?hiooo times intensity i ( * ) and toon times intensity 

I 2 ( * 5 

DATA RFDY/72H1 OOO TIMES INTENSITY II*) AND 1000 TIMES INTENSITY 

17 1*) 

j DATA CF0V/77H 1000 TIMES l INTENSITY 1 ♦ INTENSITY 2 ) 

1 

DATA OF DY/77H POLARIZATION 

t 1 

COMMON N.XN.ALP.BET.AUX1 . AUX 2. AUX3 . AUX4. AUX5 (7 ) . AUX6 (2 ) . AUX. 

1 AUKS » VX 

COMMON VI41.IST0REI901) 

5 CALL IDENT(l) 

1 FORMAT 114) 

7 F0RHATI3F10.3.5X.3I1 ) 

3 FORMAT 1 1HD»’ REFRACTIVE INDEX REIM) = ’, F7 .3 .3X. *IMl M) = »,F7.3. 

1 /» 1 H , "PARAMETER ALPHA r’,F7.2.Z) 

4 FORMAT (F8» 3) 

5 F0RMATI1H0.19X.’ AMPLITUDE FUNCTIONS’ .//. 3X .’ N’ »7X. ’REl AN I • » 
18Xt’IMIAN) , »8X*’RE(BN) ’,8X.’IMIBN) ’./) 

( . FORMAT 114. 4F 14.6) 

7 FORMAT 11 HD »9X. ’EFFICIENCY FACTOR FOR EXTINCTION OEXT =• . 
1F8.5.//.10X. ’EFFICIENCY FACTOR FOR SCATTERING OSCAT =» . 
1F8.5.//.10X. ’EFFICIENCY FACTOR FOR ABSORPTION OARS r« , 

1F8.5./) 

8 F0RMATI1HD*’SCAT ANGLE ».8X.’ INTENSITY FUNCTIONS’ »4 X. 

1 ’DEGREE OF POL ARIZATION’ ./» 1 HO »* THETA! DEG1 ’. 9X . 

1 ’ll’ .13X»*I2*.I4X, ’P’,13X.’I1*I2’./) 

9 FORMAT 1 1 H .I6.4X.2F15.6.4X.F9.6.F15.6) 

10 FORMAT 1 1 H .F7. 3. 3X .?F1 5. 6. 4X »F9. 6.F1 5. G) 

1010 FORMAT ( 1H1 .’DAT A SET NO. ’.I3.5X. »MIE SCATTERING’) 

PI=3. 1415927 
DEG-PI/180. 

READ 15 » 1 INI 


NNrl 

50 CONTINUE 

WRITE IG»10in INN 
C 50 READ!... .ETC. 

READI5»2)ALP»EMI 1) ,EMI2) .N0PT1 .N0PT2.N0PT3 
WRITE 1 6. 3) EM II ) » EM 12). ALP 
SO T015! .52) .N0PT2 

51 WRITE 16.5) 

52 CALL ZXRIEM.ALP.BET. 1) 

GO T01 13.14) .N0PT1 

♦ 13 H=37 

GO TO 53 
14 READ 15 .1 ) M 

53 DO 11 I - 1 » M 

GO T0I54.5S) .N0PT1 

* 54 XTH1 1 ) rFLO AT II— 1 ) * 5. *DFG 

GO TO 56 

55 REACTS. 41XTHII) 

XTHII)=XTHII)*DEG 

56 SIR 1 1 ) rO * 

SI II 11=0. 

S2Rll)=0. 

11 S2I I T) -0. 

Nr 1 


AEXT=D. 

ASCATrO. 
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PROGRAM LIST 3 . (Continued) 


12 XNrFLOAT ?M ) 

CALL RECJ?XJM*XJNPeXJNl« X JN1 P* X JNB » X JNBP * AUS1* AUS2.AUS3J 

AUXlrXJN 

AUX2-X JNt 

AUX3=AUS1 

AUX4=AUS2 

CALL ?EQUA?AUX5»XJNB> 

CALL ZEQUAf SUX6*AUS3> 

RBS=ALP*XJN 

R8SP-X JN+ALP«X JMP 

CALL ZHULT?8ET«XJNB*RBSB) 

CALL 7HULT?BET»XJNBP*P1) 

CALL 7A00 ? X JNB tP 1 9 RBSBP S 
RHNK J! JrRRS 

RHNK?2S = f-i<> S**N*XJNi*ALP 
RHNKPM > rPBSP 

RHNKP? 25=! -1 r ,)**N®XJNlP*ALP*?-U**N*XJNl 
CA = ? SXN*XM*1 . J * <-l 5**10/ I XN*XN-*XM5 

CB-— C A 

CALL 7XR?RRSBP*RBS*PI*11 
CALL ZXR?R8SB«RBSP»P2sl) 

CALL ZHIILT SFHsP?$P3) 

CALL 7SUB?P1 .P3.XN0MA) 

CALL 7MULT1FM.PS »P3) 

CALL ZSUB JP3»P2*XNUMB) 

CALL ZMULT ?RHNK « R8SBP® PI J 
CALL 7MULT ?RBSB« RHNKP, P2 ) 

CALL 7MULT?FM®P2»P3) 

CALL 7SUR<P1 ®P3«XOEHAJ 
CALL ZHULTfFM»PirP3> 

CALL 7SUR?P3*P2bXD£NB) 

CALL 7X1 1 XNUH A »CA*P1»1) 

CALL 70WPI *X0EMA*AN5 
CALL 7XI?XNUKB*CB*P1*1 5 
CALL 70 IV? PS ®X0E!'J8«8N) 

GO TO? 16*17? *N0PT2 

16 WRITF?S96JN»Atm>» AN?2>»BN?1 S.BNI21 

17 CXNrXN*? XN-H »S**N 

AEXN=C5?N*?AN?2?-BN?2$) 

AEXT-AEXT«-AEXN 

ABSAW— 2M00 ? AM} 

ABSBN=7M00 ? BN ? 

CXN=(XN*XN*?XN->l.S*f XN + 1 . » 5/ ?XN*XN+1 .) 

ASCN=CKN»? ABSAN*ABSAN*ABSBN* ABSBN) 

ASCAT=ASCAT*ASCN 

IF C ABS AN— .>00000015 21 *21 » 23 

21 IF (ABS8W-. D000001S22»2 2. 23 

22 IF1XN— ALPS23*23s?5 
73 GO ?OU8®i015sNOPT3 

is 00 inn 1=1 

VX=— COS? XTH C I J 5 
IFIN-l 541*41 *42 
42 AUKS=AA13KSUJ 
AUX=AAUX?TS 

41 CALL RECP?XLNP»XLNPP*AUS) 

PN 1 1 5 = XLAIP 
PNP III rXLNPP 
A AUKS? II rPN ? 1 5 
AAUX?I>=AUS 

XX?I? = XN*?XN + 1«.5*PN(I)-VX*PNPI I) 

S1NR?I 5=AM!1 >*PNP?I) +8N? 1 J*XX? 15 



PROGRAM LIST 3. (Concluded) 


ion 

101 

25 


103 


115 

117 

lie 

*J05 



102 
101 1 


EL 

1012 

104 

110 


SlNICX)rAN(21*PNPII)+BN12)*XXtIl 
S2NR(I)=BN(1 )*PNP(I)*AN( 11*XX(I) 
S2NI (T):BN(2)*PNP(I)-.AN(2)*XX( II 

S!R(I)-=SlRm*SlNR(I) 

51I(I)r:SlItII*SlNI«I) 

S2R ! II rS?R (I > *S?NR ( I ) 
S?r(I)rS?I(I)+S2NI(I) 

NrN*! 


GO TO 12 

ASCAT=ASCAT* 2 ./( ALP*ALP) 

AEXT=AEXT* 2 ./(ALP*ALP) 

AABS-AEXT— ASCAT 
WRITE! 6 ® 7 ?AEXT.ASCAT»AABS 
GO TO! 103 ® 1041 »N 0 PT 3 
WRITE (£. 8 ) 

DO 10? 1-ls.M 

xii=siRm*siR!ii+sii(ii*siim 

XI 2 =S?R(I)*S 2 R(I)»S 2 1 ( 11*5 21 III 
IEIXI?I 115 » 115»116 
Xl 2=.0000001 
WRITE ( 6 . II 7 ) I 

FORMAT 3 - 1 H .SX. 9 HBRRA Y 2 (.I 5 . 32 H) IS LESS THAN 
P 0 LAP=(XIl-XI 2 )/(XIl'*XI 2 ) 

GO T 0 (I 05 » 10 e> .N 0 PT 1 


NTHr 81—1 1*5 


X 12 =XI 1 *XI 2 

WRITE ( 6 . 9 ) N 7 H. XII. XI 2 . POLAR® XI 2 


FLOAT! NTH) 
FLOAT ( NTH) 
FLOAT (NTH) 
FLOAT ( NTH) 
XII 
XI 2 
XI 2 
POLAR 


ARRAY1 CI> 

BRRAY1 CI> 

CRRAY1 81) 

DRRAY1 (I) 

ARRAY?!!) 

BRRAY? ( I ) 

CRR A Y? ( I ) 

DRRA Y? (1 8 
GO TO 10? 

AN6=XTH(I)/0EG 

X12=XI1*XI2 

WRITE!G»10)ANG.XT1.XI2»P0LAR.X12 
SO TO 10? 

CONTINUE 

GO TO U011. 1012). N0PT1 

CONTINUE 

CALL SMXYV (0.1) 

CALL QUKLOG ! -1 .0.. 180. t.l. 1000.. 1H* 
CALL QUKLOGLO.O. .180.® .1 .1000. »1H*. 
CALL QUKLOG ! — 1 .0 .. 18 0. .. 1 . 1000 .. 1H* 


. AFOX.AFOY. 
AFOX » BFDY »— 
.AFOX.CFOY. 


C ALL SMXYV 8.0 . 0 ) 

CALL 0 UIK 3 V(- 1 . 1 H*. AFOX. DFDY.-M* DRRA Y 1 . 0 RRAY 2 ) 

CONTINUE 

NN=NN*1 

IF«NN-NI)50.50rl 10 

CALL END JOB 

STOP 

END 


OR EQUAL TO ZERO. ) 


-M.ARRAY 1 . ARRAY 2 J 

M. 8 RRAY 1 .BRRAV 2 ) 

-H.CRRAY 1 .CRRAY 2 ) 



where L 


XL, XR 
YB, YT 
ISYM 


Fixed point integer which represents the number of graphs 
per frame (1, 2, or 3). If L is negative, the frame will 
be advanced, the grid for the first graph plotted, scaled, 
labeled, and the X and Y points will be plotted. If L is 
positive, the grid for the second or third is plotted, 
scaled, and labeled, and the X and Y points are plotted. 

If L is zero, the X and Y points are plotted using the 
previous grid. For example, suppose that the user wants 
three graphs on a frame. If two graphs are desired on 
the middle grid, then in four calls the L values, in order, 
would be -3, 3, 0, 3. 

Floating point data values which represent the least and 
greatest X value respectively (for horizontal axis). 

Floating point data values which represent the least and 
greatest Y value respectively (for vertical axis). 

Fixed point integer which contains the plotting symbol 
to be used. This may be a single Hollerith character or 
an integer identifying the desired symbol. The Univac 
1108 and the SC-4020 do not have identical characters. 


FLDX, FLDY Fixed point integer arrays which contain Fieldata labels 
for the X and Y axes. FLDX and FLDY must each con- 
tain 12 words. The information to be centered must be 
left justified in the 72 character field. Remaining 
characters not used must be filled with blanks. 

NP Fixed point integer which represents the number of points 

to be plotted. If NP is negative , the plotted points will 
be connected by straight line segments. 

X, Y Floating point one dimensional arrays which contain X- 

and Y-coordinates of data to be plotted. X contains 
X-coordinates. Y contains Y-coordinates. 

The very last call of SC 4020 routines must be the card CALL ENDJOB (G) , 
or some of the last frame of the graphs might be lost. 
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The lines marked by an asterisk (* ) are for the cards that can be 
altered to change the number of points plotted. For example, if the 5 in these 
cards (e, g. , XTH(l) = FLOAT (1-1)* 5* DEG) is replaced by a 1 (as it was in 
the Atlas for a = 10-*- 100), and M is changed to 181, then points will be 
plotted at 1 degree intervals of 0 instead of 5 degree intervals. It should be 
remembered that all the 37’ s in the first dimension statement (with an asterisk) 
should be changed to 181's. The results of this program in printout for 
a = 5. 0 are found in Table 1 of the text. The SC 4020 plot for a = 5. 0 is in 
the Atlas. 
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APPENDIX B 


NUMERICAL INTEGRATION METHODS USED 
Gauss-Legendre Quadrature 

Gauss' formula for integration [5] is: 


1 n 

f f (x) dx = Yj 


w. f(x.) + R 
11 n 


R is a remainder term, negligible for the applications in this report. The 
^ th 

x.'s are the i zeros of Legendre polynomials, P^(x) , and the w.'s are the 
weighting functions determined by 


w. 

l 


2 

(1 - x.) 2 



9 


and are listed in Table 4, for n = 16. 

TABLE 4. GAUSS-LEGENDRE WEIGHTS AND ABSCISSAS 



±x. 

1 


n = 16 


w. 

1 



0.09501 

25098 

37637 

440185 

0. 18945 

06104 

55068 

496285 

0.28160 

35507 

79258 

913230 

0.18260 

34150 

44923 

588867 

0.45801 

67776 

57227 

386342 

0.16915 

65193 

95002 

538189 

0.61787 

62444 

02643 

748447 

0. 14959 

59888 

16576 

732081 

0. 75540 

44083 

55003 

033895 

0. 12462 

89712 

55533 

872052 

0. 86563 

12023 

87831 

743880 

0.09515 

85116 

82492 

784810 

0. 94457 

50230 

73232 

576078 

0. 06225 

35239 

38647 

892863 

0. 98940 

09349 

91649 

932596 

0.02715 

24594 

11754 

094852 
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For an arbitrary interval (a, b) the general formula can be modified to give: 


/ «y)dy = 2 w i f(y 1 ) + R n 

a i=l 

where, 

This integration procedure was used in the aperture and distribution 
applications. 

Gauss-Laguerre Quadrature 

The Gauss-Laguerre equation for numerical integration is of the form 

oo n 

J e" X f(x)dx = Yj w i f (V 
0 i=l 


and is useful for integrations which do not have a finite upper limit. The 
equation is used in one distribution application with a semi-infinite interval 

and a distribution according to e . 

A table of w.’s and x.*s for n = 15 is presented in Table 5. 

Trapezoidal Quadrature 

The trapezoidal quadrature method was used as a check because it is 
not as accurate as the two aforementioned methods. The extended trapezoidal 
rule (without the remainder) is 
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b 

J f(x)dx 
a 


h 


f(a) 

2 


+ f(a+ 1) + ... + f(b- i) + 


f(b) 

2 


where (a, b) is the interval, f(x) is the function, and h is the increment 
of integration [5] . 

TABLE 5. GAUSS-LAGUERRE WEIGHTS AND ABSCISSAS FOR n - 15 [5] 


x. 

1 

w. 

1 

x i 

w.e 

l 

0. 093307812017 
0.492691740302 
1.215595412071 
2. 269949526204 
3.667622721751 
5. 425336627414 
7. 565916226613 
10. 120228568019 
13. 1302.82482176 
16.654407708330 
20. 776478899449 
25. 623894226729 
31.407519169754 
38. 530683306486 
48. 026085572686 

(- 1)2.18234885940 
(- 1)3.42210177923 
(- 1)2.63027577942 
(- 1)1.26425818106 
(- 2)4.02068649210 
(- 3)8.56387780361 
(- 3)1.21243614721 
(- 4)1.11674392344 
(- 6)6.45992676202 
(- 7)2.22631690710 
(- 9)4.22743038498 
(-11)3.92189726704 
(-13)1.45651526407 
(-16)1.48302705111 
(-20)1.60059490621 

0.239578170311 
0. 560100842793 
0. 887008262919 
1.22366440215 
1. 57444872163 
1. 94475197653 
2.34150205664 
2.77404192683 
3. 25564334640 
3.80631171423 
4.45847775384 
5. 27001778443 
6.35956346973 
8. 03178763212 
11.5277721009 
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APPENDIX C 


MIE SCATTERING ATLAS 


The following graphs, produced by the SC 4020 plotter, give 

1. the scattering intensities (ij and i 2 ) 

2. the total scattering intensity (i* + i 2 ) 

3. the polarization, (i t - i 2 )/ (ij + i 2 ). 

For a = 0. 5, 0. 9, 2. 0, 3. 0, 4. 0, 5. 0, 6. 0, 7. 0, 8. 0, 9. 0, and 10. 0, the . 
points are plotted for 5-degree intervals of the scattering angle. For a = 10. 0, 
20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, and 100. 0, the points are 
plotted for intervals of 1 degree. The refractive index is m = 1.33. 

At a = 0. 5 the graphs have the form of Rayleigh Scattering. As 
a increases, the symmetry becomes more and more distorted and the complex- 
ity of the scattering diagram for large spheres is seen. 

For the first set of graphs, a = 0. 5, the meaning of each curve is 
explicitly indicated. For a = 0. 5 and 0. 9 the ordinate for the intensities 
has been multiplied by 10 3 . The points for the two components of the intensity 
are labeled with an asterisk (* ) for ij and a plus (+) for i 2 . 
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